Large Amplitude Dynamics of the Pairing Correlations in a Unitary Fermi Gas 
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A unitary Fermi gas has a surprisingly rich spectrum of large amplitude modes of the pairing 
field alone, which defies a description within a formalism involving only a reduced set of degrees 
of freedom, such as quantum hydrodynamics or a Landau-Ginzburg-like description. These modes 
are very slow, with oscillation frequencies well below the pairing gap, which makes their damping 
through quasiparticle excitations quite ineffective. In atomic traps these modes couple naturally 
with the density oscillations, and the corresponding oscillations of the atomic cloud are an example 
of a new type of collective mode in superfluid Fermi systems. They have lower frequencies than the 
compressional collective hydrodynamic oscillations, have a non-spherical momentum distribution, 
and could be excited by a quick time variation of the scattering length. 

PACS numbers: 03.75.Ss, 31.15.Ew, 67.25. dt 



While it is natural to expect the presence of hydro- 
dynamic collective modes in a unitary gas [H, the ex- 
istence of collective oscillations of large amplitude non- 
linear modes comes largely as a surprise. One reason 
is that these soliton-like modes do not emerge from a 
simplified quantum hydrodynamic or Landau-Ginzburg- 
like description of these systems. The quest for reducing 
the quantum description of a complex many-body sys- 
tem of particles to a relatively small number of degrees 
of freedom is one of the ongoing efforts in all physics 
subfields. In chemistry one would like to have an accu- 
rate description of complex molecules with many atoms 
and many more electrons in terms of a few wisely chosen 
relevant degrees of freedom (rotations, vibrations, bond 
stretching, shape, etc.). In nuclear physics major re- 
search programs are based on the assumption that many 
phenomena can be described by limiting the number of 
relevant degrees of freedom to nuclear shape and pairing 
only. Other examples are the Landau- Ginzburg descrip- 
tion of superconductors in terms of a Schrodinger-like 
description for a condensate amplitude, or the effective 
action formalism in quantum field theory which aims at 
a description of the dynamics in terms of fewer degrees 
of freedom. Most of the time it is not evident a pri- 
ori that such an approach is viable, proposed derivations 
are sometimes misleading, and many approaches rely on 
intuition and phenomenological arguments. It is partic- 
ularly important to find examples of physical systems 
where simplified descriptions are valid and the limits of 
these descriptions clear. Equally important are examples 
that show unsuspected failures of such approaches. We 
have chosen to investigate pairing dynamics in a fermion 
system, as this is relevant to a number of different prob- 
lems: nuclear collective dynamics and nuclear fission in 
particular, neutron stars, the dynamics of dilute Fermi 
gases in the unitary regime, quantum hydrodynamics of 
superfluids in general, effective action description of var- 
ious strongly interacting Fermi systems from quarks to 
nuclei to cold gases, and the sometimes invoked relation 



to a Landau-Ginzburg description of such systems. 

At first we will concentrate our attention on a uni- 
form unitary Fermi gas for a number of reasons: 1) the 
properties of both homogeneous and inhomogeneous sys- 
tems are rather well known from ab initio calculations 
0, 0|; 2) the properties of the unitary Fermi gas are 
very close to the properties of dilute neutron matter, 
which can be found in the crust of neutron stars Q ; 3) 
a very accurate description of both homogeneous and in- 
homogeneous systems is available within a Density Func- 
tional Theory (DFT) extended to describe superfluid sys- 
tems [^]; 4) the Fermi gas in the unitary regime is un- 
der intense experimental scrutiny for a number of years 
now, see Refs. 5) a wide spectrum of theoretical 

approaches of varying sophistication and accuracy have 
been applied and/or developed for this system, see recent 
review [6|; 6) the weak coupling limit has been studied 
extensively within both meanfield approximation and ex- 
act treatments 0, H, B, E3] ; 7) the experimental study of 
the aspects of the pairing dynamics discussed here ap- 
pears to be feasible in the unitary re gim e, see for ex- 
ample a somewhat related experiment [111 ]. The unitary 
regime is unlike the weak-coupling regime studied in Refs. 
0TB B EH j where no one has yet suggested a practical 
experimental realization and verification. Subsequently, 
we will discuss the case of a non-uniform unitary Fermi 
gas, where pairing gap and number density oscillations 
couple. 

To set the stage, let us review various possible frame- 
works in which one can study the dynamics of a fcrmionic 
superfluid. The oldest approach is quantum hydrody- 
namics [6|, which can be derived from conservation laws: 

n + V • \vn) = 0, (la) 
— +fi[n] + V ext \=0, (lb) 

where we dropped the arguments (r,t). Above n(r,t) 
is the number density, v(r,t) oc VS(r,t) is the velocity 
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field determined by the gradient of the phase of the con- 
densate, fi[n(r, t)] is the chemical potential, and V ext (r, t) 
is the external field in which the system might reside. 
One can derive these equations also if one assumes that 
the magnitude of the condensate remains constant and 
only its phase S(r, t) can vary. Alternatively, various au- 
thors use descriptions based on a Landau-Ginzburg in- 
spired formalisms [HI, HI , when the dynamical degrees 
of freedom are related to the "condensate wave function" 
fy(r,t). In particular these equations describe the Gold- 
stone modes arising from the broken [/(l)-symmetry in 
the condensed phase. This wave function is related to 
the expectation value of the two-fermion field operators 
ty(r,t) oc (^i(r,t)i/ji(r,t)), for which one ends up with 
a Schrodinger-like equation 
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4m 
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+U(Mr, t)|)*(r, t) + V ext (r, t)*(r, t). (2) 

Here 2m is the mass of the Cooper pair and {/(^(r, t)\) 
is a Mexican hat-like potential with a minimum at the 
ground state condensate value |^ol- One can envision 
two kinds of small amplitude oscillations in these type of 
models: 1) modes along the valley of the potential, when 
only the phase of the condensate '5 changes, which cor- 
respond to the expected Goldstone modes of the broken 
[/(l)-symmetry; 2) radial Higgs-like excitation modes, 
when the magnitude of the "condensate wave function" 
iff varies. One can derive either of these descriptions as 
a small amplitude limit of the time-dependent Hartree- 
Fock-Bogoliubov (HFB) or Bogoliubov-de Gennes (BdG) 
equations [14j . 

The HFB/BdG equations can be regarded also as an 
approximation to the appropriate time-dependent DFT 
description of such systems, namely the Time-Dependent 
Superfluid Local Density Approximation (TD-SLDA), 
which will be used here. We will show that TD-SLDA 
allows for a large number of very slow excitation modes, 
which are absent in either a quantum hydrodynamic, de- 
scribed by Eqs. |T]) or a Landau-Ginzburg description of 
these systems given by Eq. ((2J. 

The (un-regularized) SLDA energy density functional 
has been introduced and discussed in Refs. @ (see these 
references for the description of the renormalization pro- 
cedure required to eliminate the ultraviolet divergences 
of the un-renormalized theory and h = m = 1): 



£{r,t) = a-+ '- 
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where a, [3 and 7 are dimensionless parameters and 
n(r,t), r(r,t) and v(r,t) are the number, kinetic and 
anomalous densities respectively expressed through the 
usual Bogoliubov quasi-particle wave function ampli- 
tudes [ltfc(r, t),Vk(r, t)], with k labeling the quasi-particle 
states. The new element in this Letter is the non-trivial 



time-dependence of all the quasi-particle wave functions, 
which formally amounts to replacing the eigenvalues with 
time-derivatives Ek — > idt in the HFB/BdG equations. 
The time-dependent density functional theory is viewed 
in general as a reformulation of the exact quantum me- 
chanical time evolution of a many-body system when 
only single-particle properties are considered [l5j| . Af- 
ter preparing the system in its ground state, we intro- 
duce a time-dependence of 7 (which controls the magni- 
tude of the pairing field A(r,i)) on a specific schedule. 
We slowly reduce 7 in magnitude to a value j s during a 
time interval to 3> f/s Fl after which we rather abruptly 
bring it back to its value at unitarity in a time interval 
5t « 0.005/e F < l/e F , see Fig. [Hand Refs. 0,11, where 
n = k 3 F /(3ir 2 ) and e F = k 2 F /2. 

This scenario could be realized experimentally by con- 
trolling the scattering length with the magnetic field as 
a function of time, see Ref . ll| . Since by changing the 



coupling constant alone one does not induce density vari- 
ations, one might argue that the time-dependence of the 
meanfield does not a play any role in the dynamics of 
the pairing field. At unitarity one has a qualitatively dif- 
ferent scenario, since both the meanfield and the pairing 
field are of the same order of magnitude, and the mean- 
field U depends strongly on the value of the pairing field, 
see Refs. 



U(r,t) = 



/3[3ir 2 n{r,t)} 2 / 3 |A(r,t)| 2 



37n 2 / 3 (r, t) 



(4) 



We will study a range of phenomena for which the num- 
ber density n(r, t) is constant in space and time, while 
A(r, t) will be constant in space only, and this fact alone 
will lead to changes in U(r, t). In principle all couplings 
in Eq. change with the scattering length as in the ex- 
periment but since only 7 changes drastically from 
the BCS limit (k F \a\ < 1) to unitarity {l/k F a = 0), 
we neglect the changes in a and /3 for now (but rein- 
state them in a trapped system), which leads only to 
minor quantitative changes. One should remember that 
under these changes of the coupling constant (s) the num- 
ber density remains constant, and after bringing 7 back 
to its value at unitarity the total energy of the system is 
also conserved. Only during the time intervals t and St 
does one changes the energy of the system. One can ex- 
cite a large variety of oscillations, and some examples of 
the collective modes we excite in this manner are shown 
in Fig. [TJ These modes qualitatively resemble those ob- 
served in the weak coupling regime 0, IE & El • 

In panels a and b of Fig. [2] we show the instan- 
taneous single-particle occupation probabilities for the 
mode shown in panel a of Fig. [1] at times when the 
pairing field is at its smallest and and its largest val- 
ues respectively in the oscillatory regime. While the oc- 
cupation probabilities are essentially identical with their 
equilibrium values when the pairing gap is almost vanish- 
ing, the corresponding distribution at the times when the 
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FIG. 1: The panels a, b and c display response of the 
homogeneous system to an initial switching time interval 
toEp = 160, 10 and 160 and values of the gap correspond- 
ing to 7s are 7^/7 = 0.005,0.005 and 0.5 respectively, where 
7 is the coupling constant in Eq. @ and Ao ~ 0.5ef is the 
gap equilibrium value, both at unitarity. 

pairing field reaches its maximum value is drastically dif- 
ferent from its equilibrium distribution. Even though the 
occupation probabilities are so different from their equi- 
librium values, the corresponding pairing field is hardly 
different in magnitude from its equilibrium value. The 
most surprising feature is the fact that the pairing field 
oscillates around mean values different from the mini- 
mum of the "effective potential" U(\^(r, t)\) in Eq. 
namely Ao ~ 0.50^ |2J. One would naively expect that 
a simpler Landau- Ginzburg-like description as used in 
Refs. (T2I. [l3l| . would perhaps be appropriate. 

An interpretation of these modes, even in the small 
amplitude limit, as being a radial- like oscillation of the 
pairing field in a Mexican hat-like potential is equally 
invalid. In the weak coupling limit this was demon- 
strated by Volkov and Kogan [lj|, who have shown 
that the oscillations of the pairing field couple with ex- 
cited quasiparticles with energies above 2Aoo, see also 
Refs. [8|, [l0|, and that leads for large times to A(i) = 
Aoo + Asin(2A oc t + <f>)/-y/K^i. Even though the TD- 
SLDA equations have a more complex structure, the dy- 
namics of these modes is very similar in our case, see 
panel c of Fig. [TJ 

The modes displayed in panels a and b of Fig. Q] are 
truly nonlinear and their frequency depends strongly on 
the oscillation amplitude, see the panels c and d of Fig. 
[21 These are at the same time very slow modes, with 
frequencies well below the pairing gap < 2A , but 
at the same time they are truly large amplitude collec- 
tive modes, not only because of the size of the oscillation 
amplitude, but also because their excitation energy is 
equally large. In this respect these modes are somewhat 
similar to the very large amplitude oceanic waves, which 
take a long time to dissipate into heat. The results of 
the MIT experiment suggests that the damping of 
these modes, due to the decay into other modes, is much 



FIG. 2: (Color online) Panels a and b display the instanta- 
neous occupation probabilities of the mode shown in upper 
panel of Fig. [JJ corresponding to times t > when the pair- 
ing field is at its minimum and maximum values respectively 
with circles joined by a solid (blue with circles) line. With 
(red) dots we plotted the equilibrium occupation probabili- 
ties corresponding to the same instantaneous values of the 
pairing gap. In panels c and d we show the maximum and 
minimum values of the oscillating pairing field and the corre- 
sponding excitation energy as a function of the frequency of 
the Higgs-like modes, see Fig. [JJa and b. 




50 100 150 



FIG. 3: (Color online) The color bar shows the correspon- 
dence between various values of the ratio n(x, t)/n(0, 0) and 
the colors used to represent them. Here n(0, 0) = k%/(3iv 2 ) 
and £f = kp/2. The solid black line shows the correspond- 
ing rms cloud radius, see Eq. ([5}. The dashed black line 
shows the quadrupole moment of the momentum distribution 
P20 = 20{kl + k 2 z - 2kl)/(Nk%) (scaled to fit in figure). 

smaller than one might naively expect at unitarity. 

Since it would be rather difficult to study a unitary 
homogeneous system experimentally, we have considered 
the effects one might observe instead in a trapped sys- 
tem. Upon the change of the scattering length the size 
of a cloud and its central density change and this will 
induce both number density and pairing gap oscillations. 
We have considered a semi-realistic case, which we could 
simulate using our present computer resources, a homo- 
geneous system in yz-spatial dimensions (Lkp ss 131), 
which is trapped only in the third dimension in a har- 
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monic potential well V(r) = lo 2 x 2 /2 (lo/ef ~ 0.0683 and 
N = 20, 000 particles). The initial state is that of a very 
weakly interacting Fermi gas in this potential well and 
at time t — we bring quickly the scattering length to 
its unitary value and keep it constant for the rest of the 
time evolution of the system. Since the equilibrium ra- 
dius of a weakly interacting Fermi gas exceeds that of 
a unitary gas, the system tends to shrink initially and 
density oscillations of the cloud are thus excited in the x- 
direction. This is unlike the homogeneous case discussed 
above, where only oscillations of the pairing field were 
induced. Fig. [3] demonstrate a rather complex cloud dy- 
namics. The dynamics of the pairing field is very similar. 
The rms cloud size in the x-direction is described rather 
accurately as simple damped harmonic oscillations 

X(t) = ^fW)(t) =xq + Xi exp(-ijtff) cos(fttft) (5) 

with {Iff/w = 1.74 and tjh/u = 0.082 for the case illus- 
trated here and similar ratios for other cases studied by 
us. The frequency of this mode is consistently lower than 
that of the quantum hydrodynamic frequency obtained 
in the small amplitude limit from Eqs. (JTJ) for a unitary 
gas, namely SIqhd/u = 4/V5 ~ 2.31. In this respect 
this is similar to the behavior of the Higgs-like modes in 
homogeneous systems. The detailed dynamics is rather 
complicated and one can identify rather easily running 
waves, the interference of which leads to "Landau" damp- 
ing, with rjH oc vf/X (X is the system Thomas-Fermi 
radius). This is similar to waves on a surface of a wa- 
ter pool, when multiple reflections from walls lead to a 
very choppy surface, before the wave energy is converted 
into heat. This dam ping mechanism is different from 
that discussed in Ref. [17( (partially already included in 
the present approach) and likely a more efficient one as 
well in traps. We estimated the speed of these running 
waves vh/vf to be within 10% of the ground state value 
of the speed of sound c/vf — \/£/3 ~ 0.37, where £ is 
the Bertsch parameter |2j. These waves propagate with 
essentially constant speed, even though the local density 
or Fermi velocity changes quite dramatically across the 
cloud. Upon crossing two waves propagating in opposite 
directions seem to retain their original form (a soliton- 
like property), in spite of significant nonlinearities. A re- 
markable feature of this new type of excitation mode of 
Fermi systems is the intrinsic non-sphericity of the Fermi 
surface (resembling Landau's zero-sound for normal sys- 
tems), a feature absent for the hydrodynamic modes or 
the collective modes in homogeneous systems described 
above within TD-SLDA. Note that in finite systems the 
equilibrium local momentum distribution is elongated 
along the density gradient [11], thus (k 2 + k 2 — 2k 2 ) < 
initially. 

In summary, we have discussed the large amplitude 
pairing field dynamics in both homogeneous and inhomo- 
geneous unitary Fermi systems and have demonstrated 
the existence collective modes with frequencies lower 



than the hydrodynamic ones, and which in traps have 
a non-spherical oscillating momentum distribution. 
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